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ABSTRACT 

We investigate the effect of dust on the scaling properties of galaxy clusters based 
on hydrodynamic iV-body simulations of structure formation. We have simulated five 
dust models plus a radiative cooling and adiabatic models using the same initial con- 
ditions for all ru ns. The numerical imple mentation of dust was based on the analytical 
computations of lMontier fc Giardl (|2004l ). We set up dust simulations to cover different 
combinations of dust parameters that put in evidence the effects of size and abundance 
of dust grains. Comparing our radiative plus dust cooling runs to a purely radiative 
cooling simulation we find that dust has an impact on cluster scaling relations. It 
mainly affects the normalisation of the scalings (and their evolution), whereas it in- 
troduces no significant differences on their slopes. The strength of the effect critically 
depends on the dust abundance and grain size parameters as well as on the cluster 
scaling. Indeed, cooling due to dust is effective at the cluster regime and has a stronger 
effect on the "baryon driven" statistical properties of clusters such as Lx — M, Y — M, 
S — M scaling relations. Major differences, relative to the radiative cooling model, are 
as high as 25% for the Lx~M normalisation, and about 10% for the Y — M and S — AI 
normalisations at redshift zero. On the other hand, we find that dust has almost no 
impact on the "dark matter driven" T mw — M scaling relation. The effects are found 
to be dependent in equal parts on both dust abundances and grain sizes distributions 
for the scalings investigated in this paper. Higher dust abundances and smaller grain 
sizes cause larger departures from the radiative cooling (i.e. with no dust) model. 
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1 INTRODUCTION 

From the first stages of star and galaxy formation, non- 
gravitational processes drive together with gravitation the 
formation and the evolution of structures. The complex 
physics they involve rule the baryonic component within 
clusters of galaxies, and in a more general context within 
th e inte r galact ic medium (IGM hereafter - see the review 
by IVoitJ ( 20051 ) and references therein). The study of these 
processes is the key to our understanding of the formation 
and the evolution of large-scale structure of the Universe. 
Indeed, understanding how their heating and cooling abil- 
ities affect the thermodynamics of the IGM at large scales 
and high redshifts, and thus that of the intra-cluster medium 
(hereafter ICM) once the gas get accreted onto massive ha- 
los is a major question still to be answered. The continu- 
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ous accretion and the merger events through which a halo 
assembled lead to a constant interaction of the IGM gas 
with the evolving galactic component. Within denser envi- 
ronments, like clusters, feedbac k provided by AGN balance s 



Cattaneo fc Tevssierl (l2007t) 
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the gas cooling (see for ins tance J 
IConrov fc Qstrikerl (|2007l h and 
for a review). Also, from high redshift, the rate of super- 
novae drives the strength of the galactic winds and thus the 
amount of ma terial that ends eje cted within the IGM and 
the ICM (see iLoewensteinl (|2006h '). These ejecta are then 
mixed in the environment by the action of the surrounding 
gravitational potential and the dynamics of cluster galaxies 
within. 

Since long, X-ray observations have shown the abundant 
presence of heav y elements within the ICM (see for instance 
review works by (|Sarazinll 19881 : IX rnau dl2005l )h Physical pro- 
cesses like ram-pressure stripping, AGN interaction with the 
ICM, galaxy-galaxy interaction or mergers are scrutinized 
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within analytical models and numerical simulations in or- 
der to explain the presen c e of metals (see for ins t ance works 
bv (|Kapferer et al 1 l2006l : iDomainko et al.ll2006l : iMoll et al] 
[20071)) . Moreover, it is obvious that the process of tearing 
of material from galaxies leads not only to the enrichment 
of the ICM/IGM in metal, but in gas, sta rs and dust as 
well. Recent work on numerical simu lations (|Murante et al] 
|2004 120071 ; IConrov fc Ostrikerll2007l ) have stressed the role 
of hierarchical buiding of structures in enriching the ICM 
with stars in a consistent way with the observed amount 
of ICM globular clusters, and ICM light. Indeed, the overall 
light coming from stars in between cluster galaxies represent 
an important fraction of the total cluster light: for instance 
l|Krick fc Bernst ein 2007) measured 6 to 22% from a sample 
of 10 clusters. The effect of a diffuse dust component within 
the IGM/ICM, and its effect is less known. A few observa- 
tional studies with the ISO and the Spitzer satellites have 
tried witho ut frank success to detect the signature of such a 
component JStickel et al.lll998l,l2002l;jBai et alj|2006l . l2007l ). 
More successfully, ( Montier fc Giardll2005l ) have obtained a 
statistical detection, via a stacking analysis, of the overall 
infrared emission coming from clusters of galaxies. However, 
they were not able to disentangle the IR signal from dusty 
cluster galaxies from a possible ICM dust component. On 
the other hand, from a theoretical point of vi ew a few works 
have looked at the effect of dust on the ICM |Popescu et al] 
|2000| ) o r in conjunction wit h the enrichment of the ICM in 
metals |Aguirre et al.|[200ll ). However, the effect of dust on 
a I CM/IGM- type thermali zed plasma has been formalized 
by (|Montier fc Giardl I2004T I. These authors have computed 
the cooling function of dust taking into account the ener- 
getic budget for dust. They have shown the ability of dust 
to be a non negligible cooling/heating vector depending on 
the physical properties of the environment. 

Dust thus comes, within the ICM/IGM, as an added 
source of non-gravitational physics that can potentially in- 
fluence the formation and the evolution of large scale struc- 
ture in a significant way. Indeed, since redshift of z ~ 2 — 5 
during which the star formation activity reached its maxi- 
mum in the cosmic history, large amounts of dust has been 
produced and thus ejected out of t he galaxies due to vi- 
olent galactic winds into the IGM (|Springel fc Hernquistl 
I2003T I. As this material is then accreted by the forming ha- 
los, one can wonder about the impact produced by dust on 
the overall properties of clusters of galaxies once assembled 
and thermalized. In a hierarchical Universe, the population 
of clusters is self-similar, thus is expected to present well 
denned structural and scaling properties. However, to date, 
it is common knowledge that the observed properties devi- 
ate form the prediction by a purely gravitational model (see 
(|Voitll2005l : lArnaud|[200a ) for review works). It is thus im- 
portant to address the issue of the impact of dust on the 
statistical properties of structures such as clusters of galax- 
ies, the same way it is done for AGNs, supernovae, stripping 
or mergers. 

In order to tackle this question, we have put into place 
the first N-body numerical simulations of hierarchical struc- 
ture formation implementing the cooling effect of dust ac- 
cording to the dust nature and abundance. In this paper, we 
present the first results of this work focusing at the scale of 
galaxy clusters, and more specifically on their scaling prop- 
erties. The paper is organized as follows: we start by pre- 



senting the physical dust model and how it is implemented 
in the numerical simulation code. In Sec . [3] we describe 
the numerical simulations and the various runs (i.e. model) 
that have been tested. From these simulations our analysis 
concerns the galaxy cluster scale, and focus on the impact 
of the presence of dust on the scaling relation of clusters. 
In Sec. [4j we present our results on the M — T, the S — T, 
the Y — T and the Lx — T relations. The derived results are 
presented in Sec. [S] and discussed in SecOJ] 



2 THE DUST MODEL 

In our numerical simulations the implementation of the 
ph ysical effect of dust gra ins is based on the computation 
by Montier fc Giardl (|2004T ) of the dust heating/cooling func- 
tion. In this work, we decided to limit our implementation 
to the dust cooling effect only. Indeed the goal of this paper 
is to study the effect of dust at the galaxy cluster scales. 
The heating by dust grains is mainly effective at low tem- 
peratures (i.e T e < 10 5 K) and is a localised effect strongly 
dependent of the UV radiation field. Our numerical simu- 
lations (see Sec. [3] and I6.2|l do not directly implement this 
level of physics. 

Dust grains in a thermal plasma with 10 6 < T < 10 K 
are destroyed by thermal sputte r ing, w hich efficiency was 
quantified by iDraine fc Salpeter] (1 19791 . see their Eq. 44). 
The sputtering time depends on the column density and 
on the grain size. For grain sizes ranging form 0.001/jm to 



0.5/im, and an optically thin plasma (n ~ 10 atom/cm ), 
the dust lifetime spawns from 10 6 yr for small grains up 
to 10 9 yr for big grains. This lifetimes are therefore large 
enough for the cooling by dust in the IGM/ICM to be con- 
sidered. Evidently, it is also strongly linked to the injection 
rate of dust, thus to the physical mechanism that can bring 
and spread dust in the IGM/ICM. 

Our impleme ntation of the dust cool ing power is based 
on the model by (|Montier fc Giardl [2003 '). We recall bellow 
the main aspects of this model and describe the practical 
implementation within the iV-body simulations. 



2.1 The dust cooling function 

Dust grains within a thermal gas such as the ICM or the 
IGM can either be a heating or a cooling vector depend- 
ing on the physical state of the surrounding gas and on the 
radiative environment. Heating can occur via the photo- 
electric effect if the st ellar radiation field (stars and/or 
QSOs) is strong enough (|Weingartner! (|2006D and references 
therein). Indeed, the binding energies of electrons in dust 
grains are small, thus allowing electrons to be more easily 
photo-detached than in the case of a free atom or a molecule. 
On the other hand, the cooling by dust occurs through re- 
radiation in the IR of the collisional energy deposited on 
grai ns by impinging free elec trons of the ICM/IGM 0. 

iMontier fc Giardl l| 20041 ') have computed the balance of 
the heating and cooling by dust with respect to the dust 
abundance: cooling by dust dominates at high temperatures 



In the galactic medium the cooling occurs through re-radiation 
of the power absorbed in the UV and visible range. 
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in the hot IGM of virialized structures (i.e clusters of galax- 
ies), and heating by dust dominates in low temperature 
plasma under high radiation fluxes such as in the proximity 
of quasars. The details, of course, depend on the local phys- 
ical parameters such as the grain size and the gas density. 

Assuming local thermal equilibrium for the dust, the 
overall balance between heating and cooling in dust grains 
can be written as follows: 



A B (a,T d ) = H° oll (a,T e ,n e ), 



(1) 



with Hcoii being the collisional heating function of the grain 
and A the cooling function due to thermal radiation of dust. 
a is the grain size, T e and n e are respectively the electronic 
temperature and density of the medium and T d is the dust 

grain temperature. 

The heating of the dust grain was taken from iDwekl 

l|l98lf ) and can be expressed in a general way as: 



(2) 



where the values of a and j3 are dependent of the value of 
the ratio a 2/3 /T e . 

The relevant dust parameters affecting the cooling func- 
tion are the grain size and the metallicity. Indeed, the 
smaller the grains and the higher the metallicity, the higher 
is the cooling power of the dust. Thus the total cooling func- 
tion due to a population of dust grains can be expressed as 
a function of these two parameters as: 



A(a,T d ) = 



A s (a,T d ) 



dN{a, Z, V) 
dVdadZ 



dVdadZ 



(3) 



where dN(a, Z, V)/dVdadZ is the differential number of 
dust grains per size, metallicity and volume element. 

Cooling by dust happens to increase with the square 
root of the gas density, whereas the heati ng by dust is pro- 
portio nal to the density. As stressed by iMontier fc Giardl 
|2004r ) the cooling by dust is more efficient within the tem- 
perature range of 10 6 < T < 10 8 K (i.e 0.1 < kT < 10 keV), 
which is typically the IGM and ICM thermal condition s. 

We redirect the reader to IMontier k. Giardl (|2004 ) for 
a full description of the dust model, and a comprehensive 
physical analysis of the effect of dust in a optically thin 
plasma. 



2.2 The dust abundance 

The abundance of dust is a key ingredient to properly 
weight in our implementation. Observations indicate that 
dust represents only a tiny fraction of the ba ryonic mat- 
ter: Mdv, s t/Mg as « 0.01 in our Milky Way (|Dwek et all 
ll990T >, and this is possibly lower by a factor 100 to 1000 in 
the ICM: M d u St /M aas = 10~ 5 - 10~ 4 (|Popescu et al.ll2000l : 
lAguirre et al.|[2~00lh . We defined the abundance of dust as 
the ratio of the dust mass with respect to the gas mass: 



Z d = 



M„ as 



fd -5— Z dQ 

ZjQ 



(4) 



where Z is the metallicity in units of solar metallicity, 
Z d Q = 0.0075 is the solar dust abund ance, i.e the dust - 
to-gas mass ratio in the solar vincinity (|Dwek et al.lll990l ). 
and f d is the abundance of dust in the ICM in units of solar 
dust abundance. 



Dust enrichment occurs via the feedback of galaxy for- 
mation and evolution in the ICM through interaction, strip- 
ping, mergers, galactic winds and AGNs outburst. At all 
redshifts, it is linked to the SFR which drives the produc- 
tion of dust in cluster galaxies. However, in our hydrody- 
namic simulations (see Sect. |3J) the SFR is not physically 
modeled, but it is inferred by the cooling state of the gas 
particles within the simulations: gas particles below a given 
threshold of temperature and above a given threshold of 
density are considered as colisionless matter, forming stars 
and galaxies (see Sec. [3j . In order to tackle this problem, 
we choose to directly link the dust abundance to the metal 
abundance using Eq. Q. Therefore, the dust distribution in 
our simulations mimics the metal distribution. 



2.3 Implementation in the iV-body simulations 

From the equations presented in the previous sections, we 
computed the dust cooling function according to the embed- 
ding medium temperature and (global) metallicity. In simu- 
lations, once the metallicity and temperature are known, a 
and f d are the only two parameters driving the dust cool- 
ing rate (i.e A(a, Z) = A(a, f d )). In the top panel of Fig. [T] 
we present dust cooling rates (red lines) for f d = 0.1 and 
a — 10 -3 /im (model Dl, see below) at different values of 
metallicity. T he blue and black line s are t he radiative cool- 
ing rates from lSutherland fc Dopital (|l993h and the total (i.e 
radiative plus dust cooling) rate, respectively. 

Together with an adiabatic run (i.e model A) and a 
"standard" radiative run (model C - see Sect. [3] for further 
details), we ran a total of five runs implementing various 
population of grains (i.e named Dl to D5) characterized by 
their sized and dust-to-metal mass ratio: 

• We tested three types of sizes: two fixed grain sizes 
with a — 10~ 3 /an and a = 0.5 fim), respectively labeled 
small and big. The third assumes f or the IGM dust gr ains 
a distribution in sizes as defined bv lMathis et al.l (| 19771 ) for 
the galactic dust: N(a) cc a -3 ' 5 within the size interval of 
[0.001, 0.5] ^m. It is hereafter referred as the 'MRN' distri- 
bution. 

• We investigate three values of f d : 0.001, 0.01 and 0.1. 
The two extreme values roughly bracket the current the- 
oretical and observational constraints on dust abundance 
in the ICM/IGM (i.e 10.~ 5 and 10~ 3 in terms of dust-to- 
gas mass ratio) JPopescu et ah I l2000l: lAguirre et al.l l200ll: 



gas m ass ratio) (ropcs cu et aljl^ uuu: Agmrrc et a u \zv\JU 
IChelouche et al.ll2007l : lMuller et al.ll2008l : lGiard et al.ll2008l ) 



Tab.[T]lists code names and simulation details of all runs 
used in this work. In case of models Dl to D5, simulation 
cooling rates are given by the added effect of cooling due to 
dust and radiative gas cooling. Total cooling functions are 
displayed (non-coloured lines) in the bottom panel of Fig [T] 
for each of these models at Z/Zq = 0.33. As the Figure indi- 
cates, the effect of dust cooling is stronger for models with 
higher dust-to- metal mass abundance parameters, f d , and 
for smaller grain sizes (model Dl). For low values of f d the 
impact of dust cooling is significantly reduced. For exam- 
ple, in the case of model D5, the contribution of dust to the 
total cooling rate is negligible at Z/Zq — 0.33 for all tem- 
peratures. Therefore we do not expect to obtain significant 
differences between simulations with these two models. 
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Figure 1. Cooling functions implemented in the numerical sim- 
ulations. Top panel shows the dependence of dust model Dl 
(fd = 0.1 and a = 10 -3 /im) with metallicity (and tempera- 
ture) whereas the bottom panel shows different dust models at 
the same metallicity Z/Zq = 0.33 (see text). Black, blue and 
red curv es are the total cooling funct ions, radiative cooling of the 
gas from lSutherland fc Dopital l|l993h and dust cooling functions, 
respectively. 



3 NUMERICAL SIMULATIONS 
3.1 Simulation description 

Simulatio ns were carried ou t with the public code package 
Hydr a, (|Couchman et all 1 19951 ; IPearce fc Couchmanl 
Il997l ). an adapt i ve particle-particle/particle-mesh 
(AP 3 M), ij Couchmanl | 199ll ) gravity solver with a for- 
mul ation of smoothed p a rticle hydrodynamics (SPH), 
see iThacker fc Couchmanl (|200fj ). that conserves both 
entropy and energy. In simulations with cooling gas par- 
ticl es are allowed to co ol usin g the method described 
in iThomas fc Couchmanl (|l992l ) and the cooling rates 
presented in previous Section. At a given time step, gas 
particles with overdensities (relative to the critical density) 



Run 


Physics 


fd 


Grain size 


■^V steps 


A 


adiabatic (no dust) 






2569 


C 


cooling (no dust) 






2633 


Dl 


cooling with dust 


0.100 


small 


2944 


1)2 


cooling with dust 


0.100 


MRN 


2920 


D3 


cooling with dust 


0.100 


big 


2886 


Dl 


cooling with dust 


0.010 


MRN 


2698 


D5 


cooling with dust 


0.001 


MRN 


2633 



Table 1. Simulation parameters: f d , dust-to-metal mass ratios 
(see Eq. [2}, grain sizes, and number of timesteps taken by simu- 
lation runs to evolve from z=49 to z=0. Cosmological and simu- 
lation parameters were set the same in all simulation, as follows: 
Q = 0.3, n A = 0.7, n b = 0.0486, cr 8 = 0.9, h = 0.7, boxsize 
L = 100fe~ 1 Mpc, and number of baryonic and dark mater parti- 
cles, TV = 4, 096, 000. 



larger than 10 4 , and temperatures below 1.2 x 10 4 K are 
converted into collisionless baryonic matter and no longer 
participate in the gas dynamical processes. The gas metal- 
licity is assumed to be a global quantity that evolves with 
time as Z = 0.3(t/to)Z©, where Zq is the solar metallicity 
and t/to is the age of the universe in units of the current 
time. 

All simulations were generated from the same initial 
conditions snapshot, at z — 49. The initial density field 
was constructed, using N = 4, 096, 000 particles of bary- 
onic and dark matter, perturbed from a regular grid of 
fixed comoving size L = 100 /i _1 Mpc. We assumed a A- 
CDM cosmology with parameters, SI = 0.3, Qa = 0.7, 
n b = 0.0486, cr 8 = 0.9, h = 0.7. The amplitude of the matter 
power spectrum was normalized using as = 0.9. The matter 
power spectrum transfer function w as computed using the 
BBKS formula (|Bardeen et al.l|l986h. with a sh ape param- 
eter r given by the formula in ISugivamal (119951 ). With this 
choice of parameters, the dark matter and baryon particle 
masses are 2.1 x 10 10 h~ 1 M Q and 2.6 x 10 9 /i _1 M Q respec- 
tively. The gravitational softening in physical coordinates 
was 25 fe _1 kpc below z = 1 and above this redshift scaled as 
50(l + ,z)- 1 ^ 1 kpc. 

We generate a total of 7 simulation runs, listed in Ta- 
ble [1] The first two runs, which will be referred hereafter as 
'adiabatic' (or model 'A') and 'cooling' (or model 'C') simu- 
lations, do not include dust. Simulations 3 to 7 differ only on 
the dust model parameters assumed in each case, and will 
be referred to as 'dust' runs, and are labeled as 'Dl' to 'D5' 
models (see Sec. 12.31 for details on the dust models defini- 
tion). This will allow us to investigate the effects of the dust 
model parameters on our results. The last column in the 
table gives the total number of timesteps required by each 
simulation to arrive to redshift zero. For each run we stored 
a total of 78 snapshots in the redshift range < z < 23.4. 
Individual snapshots were dump at redshift intervals that 
correspond to the light travel time through the simulation 
box, ie simulation outputs stack in redshift. 
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3.2 Catalogue construction 

Cluster catalogues are generated from simulations using a 
modified version of the Sussex extra ction software devel- 
oped by Thomas and collaborators ([Thomas et al.l 1 19981; 
IPearce et al~l |2000| : iMuanwong et al. I l200ll) . Briefly, the 
cluster identification process starts with the creation of a 
minimal-spanning tree of dark matter particles which is then 
split into clumps using a maximum linking length equal to 
0.5 A^ 1 ^ 3 times the mean inter-particle separation. Here Ab 
the contrast pred icted by the spherical collap se model of a 
virialized sphere (|Eke. Navarro fc Frenklll998l ). A sphere is 
then grown around the densest dark matter particle in each 
clump until the enclosed mass verifies 



M A (< Ra) = —Rk Ap crit (z) 



(5) 



where A is a fixed overdensity contrast, p C iit(z) = 
{3H 2 /8ttG)E 2 (z) is the critical density and E(z) = 
H(z)/Ho = \J + z)' A + S1a- Cluster properties are then 
computed in a sphere of radius -R200, ie with A = 200, 
for all objects found with more than 500 particles of gas 
and dark matter. This means that our original catalogues 
are complete in mass down to 1.18 x 10 13 /i _1 M Q . For the 
study presented in this paper we have trimmed our origi- 
nal catalogues to exclude galaxy groups with masses below 
Mum = 5 x 10 13 hT 1 Mq. In this way the less massive object 
considered in the analysis is resolved with a minimum of 
2100 particles of both gas and dark mater. Our catalogues 
at z—0 have at least 60 clusters with masses above Mu m . 
This number drops to about 20 clusters at z=l. 

Cluster properties investigated in this paper are the 
mass, M, mass-weighted temperature, T mw and entropy, S 
(defined as S = feT/n _2 ' 3 ), integrated Compton param- 
eter, Y (i.e roughly the SZ signal times the square of the 
angular diameter distance to the cluster), and core excised 
(50 ft _1 kpc) X-ray bolometric luminosity, Lx. These were 
compute d in the catalogues ac cording to their usual defini- 
tions, see Ida Silva et alT(|2004 ): 



M = ^m fc , 

k 

J- mw ■ — r — 5 



2/3 



Y = 



fc B o-T (1 + x) 

m c c 2 2niH 

i 

rrii pi A bo i(Ti, Z) 



(pmu) 



(6) 
(7) 
(8) 
(9) 
(10) 



where summations with the index i are over hot (T, > 10 5 K) 
gas particles and the summation with the index k is over all 
(baryon and dark matter) particles within 7?200- Hot gas is 
assumed fully ionised. The quantities mi, Ti, n, and pi are 
the mass, temperature, number density and mass density 
of gas parti cles, respectively. Aboi is the bolometric cooling 
function in ISutherland fc Dopital (| 1993h and Z is the gas 
metallicity. Other quantities are the Boltzmann constant, 
fca, the Thomson cross-section, or, the electron mass at rest, 



m e , the speed of light c, the Hydrogen mass fraction, X = 
0.76, the gas mean molecular weight, /x, and the Hydrogen 
atom mass, ttih- 



4 SCALING RELATIONS 

In this paper we investigate the scalings of mass-weighted 
temperature, T mw , entropy, S, integrated Compton param- 
eter, Y and core excised X-ray bolometric luminosity, Lx, 
with mass, M. Taking into account Eq. © these cluster 
scaling relations can be expressed as: 



T mw = Atm (M/M ) q ™ (1 + zf™ E(z) 2/3 , 
S = Asm (Af/M ) QYT (1 + zf^ E(z)~ 2/3 , 



Y = Ayt (M/M ) aYM (1 + zf™ E(z) 



\2/5 



L X = A hM (M/M ) QLM (1 + zf™ E(z) 



\T/3 



(11) 

(12) 
(13) 
(14) 



where Mo = 10 14 /i -1 Mq and the powers of the E(z) give 
the pre dicted evolut ion, extrapolated from the self-similar 
model, (|Kaiserll 19861 ), of the scalings in each case. The quan- 
tities, A, a, and /3, are the scalings normalisation at z = 0; 
the power on the independent variable; and the departures 
from the expected self similar evolution with redshift. 

These scalings can be expressed in a condensate form, 



Vf{z) = Vo(z) {x/xo) 



(15) 



where y and x are cluster properties (e.g. T mw , M), 

yo(z)=A(l + zf, (16) 

and f(z) is some fixed power of the cosmological factor 
E{z). To determine A, a, and (5 for each s c aling we use the 
metho d described in Ida Silva et al.l l|2004 ); lAghanim et al.l 
|2008l ). To summarize, the method involves fitting the sim- 
ulated cluster populations at each redshift with Eqs. {15} 
and ()16p written in logarithmic form. First we fit the clus- 
ter distributions with a straight-line in logarithmic scale at 
all redshifts. If the logarithmic slope a remains approxi- 
mately constant (i.e. shows no systematic variations) within 
the redshift range of interest, we then set a as the best fit 
value at z = 0. Next, we repeat the fitting procedure with 
a fixed to a(z = 0) to determine the scaling normalisation 
factors yo(z). This avoids unwanted correlations between a 
and yo(z). The r.m.s. dispersion of the fit is also computed 
at each redshift according to the formula, 



(17) 



where y' = yf (see Eq. {T5}) and y[ are individual data 
points. Finally, we perform a linear fit of the normalisation 
factors with redshift in logarithmic scale, see Eq. (|16p . to 
determine the parameters A and (5. 

We note that above z = 1.5 the number of clusters in 
our catalogues decreases typically below 10, hence, we do 
not fit the scaling relations above this redshift value. 
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Figure 2. Cluster scalings at rcdshift zero for the T mw — M (top left panel), S — M, (top right panel), Y — M (bottom left panel), and 
Lx — M (bottom right panel). Displayed quantities are computed within -R20O1 the radius where the mean cluster density is 200 times 
larger than the critical density. Blue colour and triangles stand for the cooling (C) run, cyan and diamonds are for the D4 run, yellow 
and filled circles are for clusters in the D2 run, and red and crosses are for the Dl run. The lines in the embedded plots are the best-fit 
lines to the cluster distributions and the shaded areas are the fit r.m.s. dispersions for the C model, for each scaling. 



5 RESULTS 

5.1 Scaling relations at z = 

In this section we present cluster scaling relations obtained 
from simulations at redshift zero. We investigate the four 
scalings presented in Section 2] for all models under investi- 
gation. 

Figure [2] shows the T mw — M (top left panel), S mw — M, 
(top right panel), Y — M (bottom left panel), and Lx — M 
(bottom right panel) scalings, with all quantities computed 
within i?2oo- In each case, the main plot shows the cluster 
distributions for the C (triangles), D4 (diamonds), D2 (filled 
circles) and Dl (crosses) simulations, whereas the embed- 
ded plot presents the power law best fit lines (solid, triple 
dot-dashed, dashed and dot-dashed for C, D4, D2 and Dl 



models, respectively) obtained in each case, colour coded 
in the same way as the cluster distributions. Here we have 
chosen to display dust models that allow us to assess the ef- 
fect of dust parameters individually. For example, the dust 
models in runs D4 and D2 only differ by the dust-to-metal 
mass ratio parameter, whereas models D2 and Dl have dif- 
ferent grain sizes but the same fa. The shaded gray areas in 
the embedded plots give the r.m.s. dispersion of the fit for 
cooling (C) model. The dispersions obtained for the other 
models have similar amplitudes to the C case. The scalings 
of entropy and X-ray luminosity with mass show larger dis- 
persions because they are more sensitive to the gas physical 
properties (density and temperature) in the inner parts of 
clusters than the mass-weighted temperature and Y versus 
mass relations which are tightly correlated with mass. 
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An inspection of Fig. [2] allows us to conclude that the 
cluster scalings laws studied here are sensitive to the under- 
laying dust model, and in particular to models where the 
dust cooling is stronger (model Dl and D2). The differences 
are more evident in the S — M and Lx — M scalings, but are 
also visible, to a lower extent, in the T mw — M and Y — M 
relations. Generally, the inclusion of dust tends to increase 
temperature and entropy because the additional cooling in- 
creases the formation of collisionless (star forming) material 
leaving the remaining particles in the gas phase with higher 
mean temperatures and entropies. The decrease of Y and 
X-ray luminosities reflects the effect of lowering the hot gas 
fraction and density due to dust cooling. These effects dom- 
inate over the effect of increasing the temperature. 

In fact a closer inspection of Fig. [5] indicates that dif- 
ferences for the same cluster in different models (note that 
all simulations have the same initial conditions so a cluster- 
to-cluster comparison can be made), reflect the differences 
of intensity between cooling functions presented in Figure [T] 
For example, the differences between models D4 and C are 
clearly small as one would expect from the small differences 
between cooling functions displayed in the bottom panel of 
Figure [T] Another interesting example is that an increase of 
one order of magnitude in fd from D4 to D2 seems to cause 
a stronger impact in the properties of the most massive clus- 
ters than the differences arising from changing the dust grain 
sizes from D2 to Dl. Again this reflects the differences be- 
tween cooling functions, which in the latter case are smaller 
at higher temperatures (see bottom panel of Figure [TJ. 

A way of quantifying the effect of dust, is to look at the 
best fit slope, a, and normalisation, log A, parameters of 
these scalings which are presented in Table [2] for all cooling 
models considered in this paper. We find that fitting param- 
eters are quite similar for models C, D5, and D4 whereas 
models with high dust abundances provide the strongest 
variations of the fitting parameters, particularly for the nor- 
malisations. In several cases, differences are larger than the 
(statistical) best-fit errors, particularly for the Dl and D2 
models. We also investigated scalings at redshift zero for 
the A (adiabatic) model and found they were consistent 
with self-similar predictions. As expected, the results ob- 
tained for the adiabatic and cooling models are in very 
good agreement with the findings of (|da Silva et al.l [20041 ; 
lAghanim et al.|[2008t ) which use similar simulation parame- 
ters and cosmology. 

5.2 Evolution of the scaling relations 

We now turn to the discussion of the evolution of the clus- 
ter scaling laws in our simulations. Here we apply the fit 
to a power law procedure described in Section [3] to derive 
the logarithmic slope, j3, of our fitting functions, Eqs. (|lip ~ 
<|14[) . As mentioned earlier, this quantity measures evolution 
departures relative to the self-similar expectations for each 
scaling. 

In Figs. [3] |4] [5] and [6] we plot the redshift dependence 
of the power law slopes, a, (top panels), and normalisations, 
log 2/0, (middle panels) for our T mw — M, S - M, Y - M, 
and Lx — M scalings, respectively. The bottom panels show 
straight lines best fits, up to z=l, to the data points in the 
middle panels of each Figure. The slopes of these lines are 
the (3 parameters in Eqs. (|11[1 - (|14[ ). We decided not to in- 
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Figure 3. Evolution of the slope (top panel), normalisation (mid- 
dle panel), and normalisation best fit lines (bottom panel) of the 
Tinw — M cluster scaling relation for the C (triangles, solid line), 
D4 (diamonds, triple-dot-dashcd line), D3 (squares, short-dashed 
line), D2 (circles, dashed line) and Dl (crosses, dot-dashed line) 
simulation models. Colour bands are best fit errors to the clus- 
ter distributions at each redshift. The shaded area in the bottom 
panel is the rms dispersion of the normalisation fit for the cooling 
model, (see text for details) 



elude data points above z=l in the computation of j3 be- 
cause cluster numbers drop rapidly (below 20) which, in 
some cases, causes large oscillations in the computed nor- 
malizations. Moreover in the case of the Lx — M relation, 
the evolution of 2/0(2) with redshift appears to deviate from 
a straight line above z ~ 1. In Table[5]we provide a complete 
list of the log A, [3 and a fitting parameters and associated 
statistical errors for all scalings and cooling models inves- 
tigated in this paper. The displayed values are valid in the 
redshift range < z < 1. In the top and middle panels the 
coloured bands correspond to the ±1ct envelope of the best 
fit errors obtained at each redshift for a, and log 2/0 ■ The 
shaded area in the bottom panels are r.m.s. fit dispersions 
of the normalisations, log 2/0, computed for the cooling model 
using Eq. (flT)) . 

Results from different simulation runs are coded in the 
following way: triangles and solid lines stand for the cooling 
model, diamonds and triple-dot-dashed lines represent the 
D4 model, squares and short-dashed lines are for D3 model, 
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Figure 4. Evolution of the slope (top panel), normalisation (mid- 
dle panel), and normalisation best fit lines (bottom panel) of the 
S — M cluster scaling relation for the C, D4, D3, D2, Dl sim- 
ulation models. Symbols, lines and colours are the same as in 
Fig [3] 



Figure 5. Evolution of the slope (top panel), normalisation (mid- 
dle panel), and normalisation best fit lines (bottom panel) of the 
Y — M cluster scaling relation for the C, D4, D3, D2, Dl sim- 
ulation models. Symbols, lines and colours are the same as in 
Fig! 



circles and dashed lines for the D2 models and crosses and 
dot-dashed lines are for the Dl model. Here we have omitted 
the D5 model for clarity. It provides the same fitting results 
as the cooling model. This confirms expectations and the 
comments made in the last paragraph of Section I2.3I 

The top panels of these Figures allow us to conclude 
that the a slopes of our scalings are fairly insensitive to dust 
cooling. These also show no evidence of systematic variations 
with redshift for all scalings, which is an important require- 
ment when fitting the cluster distributions with power-laws 
of the form Eqs. (|1 1 fl - (| 14[l . The redshift independence of the 
slopes with the dust model confirms our findings at redshift 
zero. The scatter (non-systematic "oscillations" ) at high red- 
shift is caused by the decrease of the number of clusters with 
Mii m ) 5 x 10 13 /i _1 Mq, the sample selection used for all 
fits. 

The main effect of cooling by dust is reflected in the 
changes it produces in the normalisations of the cluster scal- 
ing laws. Again, the the impact of dust is different depend- 
ing on the scaling under consideration. For the T mw — M 
scaling in Fig. Owe see a sytematic variation with the dust 
model (ordered in the following way: C, D4, D3, D2, Dl), 
but differences between models are within the errors and fit 



dispersions of each others. For the evolution of the normal- 
isations of the S — M, Y — M, and Lx — M scalings (see 
Figs. 3] and[U} we conclude that the inclusion of dust cool- 
ing causes significant departures from the standard radiative 
cooling model depending on the dust model parameters. For 
example, this is clear from the non-overlapping errors and fit 
dispersions of the normalisations for the D2 and Dl models. 
For all scalings, the relative strength of the effect of dust fol- 
lows the relative intensity of the cooling functions presented 
in Section \2. 31 This orders the models in the following way: 
C, D4, D3, D2, Dl with increasing normalisations for the 
Tmw — M and S — M scalings and decreasing normalisations 
for the Y — M and Lx — M relations. 

We end this section by noting that we find positive evo- 
lution (relative to the expected self-similar evolution, i.e. for 
a given x in Eq {15} the property yf is higher at higher red- 
shifts) for the Y — M and L x -M (models Dl and D2 only) 
relations whereas the T mw — M, and S — M relations show 
negative evolutions relative to the self-similar model. This 
is in line with the findings from simulations with radiative 
cooling of similar size and cosmology, see eg (|da Silva et alj 
|2004 lAghanim et"aT]|2008l ) . 
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Figure 6. Evolution of the slope (top panel), normalisation (mid- 
dle panel), and normalisation best fit lines (bottom panel) of the 
Lx — M cluster scaling relation for the C, D4, D3, D2, Dl sim- 
ulation models. Symbols, lines and colours are the same as in 
Fig [3] 



6 DISCUSSION 

6.1 Efficiency of the dust cooling 

In ag reement with the cooling functions of (|Montier fc Giardl 
120041 ). the dust cooling is most effective in the cluster tem- 
perature regime. The relative importance of the dust cool- 
ing with respect to the gas radiative cooling is strongly de- 
pendent on the dust abundances and the intrinsic physical 
properties of the dust. This is clearly shown in our scaling 
relations results: 

• the T mw — M relation is almost unchanged when adding 
dust cooling to the radiative gas cooling (see Fig. |3J). Our 
results show that the (mass weighted) temperature-mass re- 
lation within -R200, is essentially driven by the gravitational 
heating of the gas (due to its infall on the cluster poten- 
tial well), and that the physics of baryons (at least for the 
physics implemented in the simulations presented in this pa- 
per) play very little role in the outer parts of halos which 
dominate the estimation of the mass-weighted temperature 
and integrated mass. Since gas cooling tends to disturb the 
dark matter distributi on at the centre of clusters in high 
resolution simulations (jGnedin et al.l [20041 ). the cooling by 
dust may amplify this effect, and thus modify scaling rela- 



tions like the T mw — M. In the case of observationally derived 
quantities, scaling laws will be drawn from overall quantities 
that will proceed from mixed-projected information over a 
wide range of radii. If a gradient exist in the dust effect 
towards the cluster centre, an "overall" temperature might 
bear the signature of the structural effect of dust. Anyway 
this quantification is beyond the scope of this paper and 
will be investigated in a forthcoming paper. There is also no 
significant effect between the different dust models and the 
radiative case on the evolution of the slope and normalisa- 
tion of the T mw — M relation. 

• On the other hand the other three scaling laws are 
deeply related to the clusters baryonic component. The clear 
effect on the S — M, Y — M, and Lx — M relations illustrates 
this fact (see Pigs. [4j [5j and [6}. We found that the slopes of 
these scalings remain fairly insensitive to dust, whereas nor- 
malisations show significant changes depending on the dust 
parameters. Relative changes in the normalisations at red- 
shift zero and Mo = W 14 h~ 1 M@ can be as high as 25% for 
L x - M and 10% for the S — M, Y — M relations for the Dl 
model. Models with lower dust abundances and MRN grain 
size distributions present smaller but systematic variations 
relative to the C model. As any other cooling process, the 
cooling due to dust tends to lower the normalisations of the 
Y — M, and Lx — M scalings due to the decrease of the hot 
gas fractions and densities which dominate the increase of 
temperature. The increase of normalisations for the S — M 
and T mw — M relations with added dust cooling is also in 
line with expectations because cooling converts cold, dense, 
gas into collisionless star forming material, which raises the 
mean temperature and entropy of the remaining gas. 

• Our simulations allow us to quantify the relative impact 
of the dust parameters on the investigated cluster scalings 
(see Figs. [3] to [5] and Table [5}. From one model to another 
one can identify two clear effects due to dust: (i) it shows the 
expected effect of the dust abundance, which from models 
D4 to D2 raise by a factor of 10, producing a change of 
normalisations relative to the purely radiative case (model 
C), from almost zero percent contribution to about 14%, 5% 
and 6% contributions for the L x - M, Y - M and S - M 
relations, respectively, (ii) Even more striking is the effect of 
the intrinsic dust grain physical properties. The variation of 
normalisations relative to the C model, change from a zero 
percent level for model D4 to about 25% (L x — M) and 10% 
(Y — M, and S — M) for the model Dl (ie the relative change 
from models D2 to Dl is about 13% and 5%, respectively). 
All these percentages were calculated using normalisations 
at redshift zero and M = 10 14 h _1 M Q . Therefore the size 
of the grains comes to be an equally important parameter 
varying the efficiency of the dust cooling. The smaller the 
grain, the stronger the cooling. 

• From Figs. [3] [4] [5] and [6] one finds that differences be- 
tween normalisations become progressively important with 
decreasing redshift. This confirms expectations because 
metallicity was modelled in simulations as a linearly increas- 
ing function of time. Although our implementation of metal- 
licity should only be regarded as a first order approximation 
to the modelling of more complex physical processes (act- 
ing on scales below the resolution scale of the present set 
simulations), it would be interesting to investigate whether 
a similar effect remains (ie the effects of dust become pro- 
gressively important at low redshift) when such processes 



10 A. da Silva et al. 



are taken into account throughout the formation history of 
galaxy clusters (see discussion below). 



6.2 Limitation of the dust implementation 

In order to implement the presence of dust in the numerical 
simulations, we chose a "zero order approach": we directly 
correlated the presence of dust with the presence of metals 
under the assumption that there is no segregation in the na- 
ture of the material withdrawn from galaxies and injected 
in the IGM/ICM (metals, gas, stars or dust). However, this 
assumption suffers from limitations linked to the dust life- 
time. Indeed, dust grains strongly suffer of sputtering and 
within their lifetime they are depleted from metals which, 
contrary to dust grains, are not destroyed i.e. remain in the 
ICM/IGM. Therefore our whole analysis is to be considered 
within the framework of this assumption, and is to be under- 
stood as a basic implementation of the effects of dust with 
the objective of assessing whether dust has a significant im- 
pact on large scale structure formation, and consequently to 
quantify these effects at first order. 

Moreover, our implementation is also ad hoc. Indeed, 
beside the cooling function of dust, our implementation is 
not a physical implementation. We did not deal stricto senso 
with the physics of the dust creation and dust destruction 
processes. This would be a step further, and is yet beyond 
the scope of this paper as mention ed above. However, mak - 
ing use of the cooling function by l|Montier fc Giardll2004 ). 
we have performed a fully self-consistent implementation of 
the effect of dust as a cooling vector of the ICM/IGM. In- 
deed, on the basis of the cooling function, the implemen- 
tation encapsulates the major physical processes to which 
dust is subjected and acts as a non-gravitational process at 
the scale of the ICM and the IGM. 

As already mentioned, we directly correlated the abun- 
dance of dust with metallicity, thus to the metallicity evo- 
lution, which chosen evolution law is quite drastic: Z = 
0.3(t/to)ZQ. Indeed, if the metallicity at z = is normal- 
ized to the value of 0.3/70, it is lowered to ~ 0.2 at z = 0.5 
and ~ 0.1 at z = 1. However, other numerical works based 
on simulations including physical implementation of metal 
enrichment processes but without dust agree well with ob- 
servational constraints (mainly provided by X-ray observa- 
tions of the Fe-K line) which indicate hi gh metalicity val- 
ues, Z ~ 0.3Z a , up to redshifts above 1.0 (|Cora et alfoost 
iBorgani et alj|2008l ). This shows that, as for the stellar com- 
ponent which is already in place in galaxies when clusters 
form, the metal enrichment of the ICM/IGM has occured 
through the feedback of galaxy formation and evolution, and 
therefore it de facto strongly enriched the IGM/ICM bellow 
z = 1. It also might give hints that the high metallicity of 
clusters could be correlated to the dust enrichment of the 
IGM/ICM. Indeed, the amount of gaseous iron in galaxies 
such as the Milky Way is ~ 0.01Z©. An early enrichment of 
dust in the IGM and/or the ICM, which once sputtered will 
provide metals, could explain part of the iron abundances 
found in the ICM at low redshifts. This hypothesis seems 
to be consolidated by the few works that have investigated 
dust as a source for metals in the material stripped from 
galax ies via dynamical re moval within already formed clus- 
ters (jAguirre et alj|200ll) or via an early IGM enrichment 



at high redshift during the pea k of star formation around 
z = 3 (|Bianchi fc Ferrara] 120051 ). The latter work stressed 
that only big grains (a > 0.1/im) can be transported on a 
few lOOkpc physical scale, however leading to a very inho- 
mogenous spatial enrichment in metals once the dust grains 
are sputtered. For all these reasons, by underestimating the 
metallicity at high redshifts, we might have underestimated 
the amount of dust injected in the ICM at high redshift, 
and thus the efficiency of dust cooling when integrated from 
early epoch down to redshift zero 



7 CONCLUSION 

In this work, we have presented the first simulations of struc- 
ture formation investigating the effect of dust cooling on the 
properties of the intra-cluster medium. We have compared 
simulations with radiative plus dust cooling with respect to 
a purely radiative cooling simulation. We have shown that: 

• The cooling due to dust is effective at the cluster regime 
and has a significant effect on the "baryon driven" statistical 
properties of cluster such as Lx — M, Y — M,S — M scal- 
ing relations. As an added non-gravitational cooling process 
dust changes the normalisation of these laws by a factor up 
to 27% for the L x -M relation, and up to 10% for the Y — M 
and S — M relations. On the contrary, dust has almost no 
effect on a "dark matter driven" scaling relation such as the 
T mw — M relation. 

• The inclusion of cooling by dust does not change signif- 
icantly the slopes of the cluster scaling laws investigated in 
this paper. They compare with the results obtained in the 
radiative cooling simulation model. 

• Through the implementation of our different dust mod- 
els, we have demonstrated that the dust cooling effect at the 
scale of clusters depends strongly on the dust abundance 
in the ICM, but also on a similar proportion on the size 
distribution of dust grains. Therefore the dust efficiency is 
strongly dependent on the nature of the stripped and ejected 
galactic material, as well as the history of these injection and 
destruction processes along the cluster history. Indeed the 
early enrichment of dust might provide an already modi- 
fied thermodynamical setup for the "to-be-accreted" gas at 
lower redshifts. 

The setup of our simulations and the limitation of our 
dust implementation can be considered at a "zero order" 
test with which we demonstrated the active effect of dust 
on structure formation and especially at the cluster scale. 
In order to go one step further, a perspective of this work 
will be needed to couple the radiative cooling function of 
dust with a physical and dynamical implementation of the 
creation and destruction processes of dust in the IGM/ICM. 
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Table 2. Best fir values of the parameters a, log A and j3 as well as their respective la errors. These values are valid within the rcdshift 
range < z < 1. 



Model C Model D5 Model D4 Model D3 Model D2 Model Dl 



T mw - M 

«tm 0.61 ±0.02 0.61 ±0.02 0.61 ± 0.02 0.62 ± 0.02 0.63 ± 0.02 0.63 ± 0.02 

logA TM 0.195 ±0.002 0.195 ±0.003 0.196 ± 0.002 0.197 ±0.003 0.201 ± 0.002 0.204 ±0.002 

/3 T m -0.14 ±0.01 -0.14 ±0.01 -0.14 ±0.01 -0.15 ±0.01 -0.16 ±0.01 -0.16 ±0.01 



S-M 

a SM 0.55 ±0.03 0.54 ± 0.03 0.54 ± 0.03 0.56 ±0.03 0.55 ± 0.02 0.54 ±0.02 

log Asm 2.443 ±0.002 2.444 ± 0.002 2.445 ± 0.002 2.451 ± 0.002 2.468 ± 0.002 2.488 ± 0.02 

/3 SM -0.33 ±0.01 -0.34 ±0.01 -0.34 ± 0.01 -0.36 ± 0.01 -0.40 ± 0.01 -0.42 ± 0.01 



Y — M 

o YM 1.74 ±0.03 1.72 ±0.03 1.73 ± 0.03 1.72 ±0.02 1.74 ±0.02 1.76 ± 0.02 

logA YM -5.909 ±0.002 -5.907 ±0.002 -5.910 ± 0.002 -5.914 ± 0.002 -5.933 ± 0.002 -5.957 ±0.002 

/3 Y m 0.12 ±0.01 0.11 ±0.01 0.13 ±0.02 0.13 ±0.01 0.17 ±0.01 0.21 ± 0.01 



L x -M 

a LM 1.69 ±0.07 1.68 ±0.07 1.65 ± 0.07 1.61 ±0.08 1.67 ±0.05 1.67 ±0.05 

log Alm 3.330 ±0.006 3.334 ± 0.006 3.333 ± 0.005 3.323 ± 0.005 3.265 ± 0.005 3.207 ±0.004 

/3 LM 0.01 ±0.03 -0.02 ±0.03 -0.02 ± 0.03 0.02 ± 0.03 0.18 ± 0.03 0.23 ± 0.03 



